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We study networks of noisy phase oscillators whose nodes are characterized by a random degree 
counting the number of its connections. Both these degrees and the natural frequencies of the oscil- 
lators are distributed according to a given probability density. Replacing the randomly connected 
network by an all-to-all coupled network with weighted edges, allows us to formulate the dynamics 
of a single oscillator coupled to the mean field and to derive the corresponding Fokker-Planck equa- 
tion. From the latter we calculate the critical coupling strength for the onset of synchronization as 
a function of the noise intensity, the frequency distribution and the first two moments of the degree 
■ distribution. Our approach is applied to a dense small-world network model, for which we calculate 

the degree distribution. Numerical simulations prove the validity of the made replacement. We also 
test the applicability to more sparsely connected networks and formulate homogeneity and absence 
of correlations in the degree distribution as limiting factors of our approach. 

PACS numbers: 05.40.-a, 05.45.Xt, 87.10.Ca, 87.19. Ij 



I. INTRODUCTION 

Synchronization phenomena have been observed in many different fields such as the chaotic 
intensity fluctuations of coupled lasers flashing swarms of flreflies 0, spiking dynamics of 
neurons pacemaker cells in the heart 0], and the menstrual cycles of women living together 

(sl-fioj , to mention only a few examples. Without doubts, the phenomenon of synchronization is a 
central mechanism in physics, chemistry, biology, and medicine as reported impressively in various 
monographs pdl - fT^ . 

The key aspect of our analysis is the phase description of the dynamics tlSi]- With regard to 
biological oscillators, the idea of the phase description goes back to Winfree (1967) who motivated 
it by the observation that the oscillators should be weakly coupled so that no oscillator is ever 
perturbed far away from its limit cycle and therefore amplitudes could be considered as flxed [l6| . 

Kuramoto (1975, 1984) simplifled the Winfree model by certain methods of perturbation and 
averaging [13, [3 ■ The Kuramoto model has been shown to be a very successful approach to the 
problem of synchronization 19]. By virtue of its simplicity the Kuramoto model is analytically 
tractable. We take advantage of the simplicity to the greatest extent, but at the same time we 
are determined to make the model more realistic, in particular with regard to neural networks. 
Therefore, we add two features: (i) we consider complex networks instead of all-to-all connectivity 
and (ii) we study noise in the phase model. The flrst feature takes account of the fact that many 
real- world networks, such as neural systems, are complex networks "par excellence" [20| . while 
the second feature incorporates experimentally observed stochastic processes, such as the random 
synaptic inputs from other neurons [isj . 
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For an ensemble of Kuramoto oscillators, both features have been considered separately in the 
literature: complex networks in [2l| - [23l a nd noise in (23 - [26| . For an overview on nonlinear dynamics 
in complex networks we recommend [27H30| and for an overview on the Kuramoto model we further 
recommend [31] 



Referring to 



33|. 

22l |. we apply an approximation technique, in which we focus on the degrees of 



the nodes, i.e. the numbers of connections leading to and emanating from a node. The degree 
distribution could be given "a priori" and then defines the network or it is determined by counting 
the frequencies of the edges for every node of the large network under consideration. 

As a result of the approximation, the complex network is replaced by an all-to-all coupled one 
with weighted edges, which mimics the original complex structure. This procedure originates a 
mean-field like description that enables us to derive both the effective Langevin equation for a 
single oscillator inside the random network and the corresponding nonlinear Fokker-Planck equation 
(FPE). Moreover, the latter gives us the possibility to calculate the critical coupling strength that 
marks the onset of synchronization in the network. 

The work is organized as follows. In section [TT] we present our extended Kuramoto model and in 
section mil we explain the approximation technique. In section [TVl we derive the nonlinear FPE and 
in |V] the critical coupling strength. This result is compared with numerical simulations on small- 
world networks (section |Vl| , and the applicability to other networks is discussed (section |V11[) . In 
the last section IVIIII we give concluding remarks. 



II. CONCEPTION OF OUR MODEL 

We refer to the original Kuramoto model [13, [3 in the following way: 

N 

Mt) = ^« + ^ H sin ('/'j(i) - Mt)) + Ut), I = 1, . . . , iV , (1) 

with the number of oscillators iV, the phase 4>i{t) at time t and natural frequency cui of oscillator 
i. The coupling strength is denoted by >c and Aij are the elements of the adjacency matrix. If the 
nodes i and j are not connected, then Aij = 0, otherwise Aij = 1. Self-coupling is excluded in 
weighted networks. The adjacency matrix further yields the individual degrees: 

N 

fc, z = l,...,7V. (2) 

The natural frequencies and individual degrees are distributed according to a joint probability 
density P(w, fc). 

Remember that the indegree or the outdegree of node i are defined as the number of edges 



pointing to or emanating from node i, respectively [34| . Since we consider undirected networks, the 
adjacency matrix is symmetric and the indegrees are equal to the outdegrees. Hence, we refer only 
to the degree ki of node i. 

The division by N in the coupling term of Eq. ([T]) is not always the correct normalization for 
complex networks, because it might not lead necessarily to an intensive coupling term. That is 
the case, if the adjacency matrix does not scale with the system size. Then one has to take the 
maximum degree instead [s^]. 
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The functions S,i(t), i — I, . . . , N, stand for sources of independent white noise processes that act 
on the natural frequencies as a stochastic force. They satisfy 

Hence, a single nonnegative parameter D scales the noise. It is its intensity and we assume that it 
is independent of the system size. The angular brackets denote an average over different realizations 
of the noise. 



III. APPROXIMATION TECHNIQUE 



The idea is to adopt a combinatorial point of view in order to mimic the complex network 
by a weighted fully connected network with a new adjacency matrix Aij which should resemble 
the original structure defined by the given set of degrees fc^, i = 1,...,N. We require that the 
approximating weights Aij conserve the degrees of the original network (cf. Eq. (HI)), i.e. 

N N 

fti = ^ Aij and kj = ^ Aij, (4) 
j=i i=i 

where we assume again an undirected network. Eq. Q also means that the new matrix scales in 
the same way as the old one, if the system size is changed. 

In order to construct the approximating weight between nodes i and j, we treat the original 
complex network as a random network and assume that the coupling strength is proportional to 
the probability that a node with degree ki couples to a node with degree kj. If the degrees are 
uncorrelated [13, [11] , then Aij is equal to the degree hi times the degree kj normalized by the total 
number of all degrees in the network: 

Obviously, it defines a symmetric matrix since the same expression is found using the same argu- 
ments for the weight of the edge connecting the j-th with the i-th node. It is easy to see that 
the new matrix conserves the local degree structure as required in Eq. (|3|) and the scaling of the 
adjacency matrix. We illustrate our approximation in Fig. [1] 
Inserting Eq. (O into Eq. ^ yields 

>c k ^ 

</.,(0 + sin (0,(t)- 0,(0) i = (6) 

This is the approximative form of our model, which appears to be analytically tractable. 
Let us make the following definition (cf. Ref. [2^): 



Ef^,fc,-e-^^ W 
k 



^(t)e^^(*' . (7) 
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FIG. 1: (Color online) The effect of our approximation on a ring network of eight symmetrically coupled 
oscillators. On the left-hand side the original complex network is shown and on the right-hand side the 
approximate network is shown, where the thickness of the edges is chosen approximately proportional to 
the coupling strength. For reasons of clarity, self-coupling is not visualized. 

The quantity r{t) is an order parameter, because for a population of A'^ — >■ oo completely 
asynchronous oscillators, r{t — > oo) vanishes, while the onset of synchronization is marked by 
r{t — >■ oo) > (the completely phase synchronized state corresponds to r{t — >■ oo) = 1). 

By multiplying Eq. ([7]) by e"*"^'^*^ and by considering only the imaginary parts, we can rewrite 
Eq. (O as an effective one-oscillator description, where the common time-dependent phase Q{t) 
and amplitude r{t) are averaged over all the nodes according to Eq. ([7]): 

Mt) + r{t)x^ sin(e(i) - (/.,(t)) + i,{t), z = 1, . . . , iV. (8) 

In other words, our approximation is equivalent to a mean- field approximation with the mean field 
amplitude r{t) and phase 8(t) defined in Eq. ([7]). All oscillators are statistically identical and 
differ by u>i and ki only. They are coupled to the mean field with a characteristic strength that 
is proportional to the individual degree ki, which can be seen as a property of node i. Such an 
assignment is similar to a characterization of subpopulations through coupling strengths [36l | . 

IV. DERIVATION OF THE NONLINEAR FOKKER-PLANCK EQUATION 

The solution of the system of TV coupled Langevin equations in Eq. © is a Markov process. It 
can also be formulated by the help of the transition probability V (0, t\cff' ,t^\ lj, fe) that describes 
the evolution of the phases from time t° to time t > . Therein, respectively the vector 4> = 
(01, . . . , (t>M) contains the phases of the N oscillators at time t and 0" = (0°, . . . , 0^) at initial 
time . The vectors uj = (wi, . . . ,uin) and k = (/ci, . . . , k^) contain the natural frequencies of the 
oscillators and the degrees of all the nodes in the network, respectively. The individual frequencies 
uj and degrees k are initially chosen at time t — and they stay fixed during the whole evolution. 
We also assume that the initial phases (j)^ are identically and independently distributed. 

The transition probability V (0, t|</)°, a;, fc) for N oscillators satisfies a linear Fokker-Planck 
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equation (FPE) 

and is subject to the initial condition 

7'(0,iO|0O,i";a;,fe) =<5^(0-0"). (10) 

Conditioned and joint probabilities satisfy 

V (0, a;, k) V (0°, t°; u,,k)^V (0, t; lj, k\4P , i°) (0°, t°) . (11) 

As already stated in section |TT1 the natural frequencies and degrees are distributed according to the 
joint probability density P{uj,k). The initial phases of the oscillators are further supposed to be 
chosen independently of a; and k. Therefore, we have the joint distribution of initial phases, edges 
and frequencies at initial time t*^ factorizing as 

N 

V (0", t°; a;, fc) = n ^ (^^ h) V t°) . (12) 

i=l 

The joint distribution at a later time t becomes independent of the initial distribution by integrating 
Eq. (fTTj) over the initial phases: 

JV 

V (0, t- u;, fe) - / d0° . . . d0^ V (0, t|0", i"; a;, fc) [] P (c.„ fc,) T' (0^, t°) . (13) 

i=i 

Again, we aim at reducing the description to an effective one-oscillator picture. As does the 
transition probability, the distribution Eq. obeys the FPE A reduced joint distribution 
Pn with n < N for n phases with frequencies wi . . . a;„ and degrees fci , . . . fc„ is defined in the usual 
way [131: 

/9„ ((/)i, . . . , 0„,i; wi, . . . ,a;„, fci, . . . , fc„) 

:= J d(j)n+i . ■ -Mn J dujn+i ■ ■ -dujN J dkn+i ■ . -dkN V {4>,t;u:,k) . ^ ^ 

For the important case of a single oscillator n ~ 1 imbedded in the network, we integrate Eq. ^ 
over the remaining — 1 phases {4>2t ■ ■ frequencies {u)2, . . . , wjv) and degrees (fc2, . . . , /cat): 



dpi _ 5 , . „ , r^d^Pi _ 

(15) 



dt - 90^^^^+ ^90? 



d 



K ki {N ~l) 



■2 / dw2 / dfc2Sin(02 - 0l)fc2P2 (01,</'2,i;Wl,W2,fci,/c2) 



where we use identity of all the oscillators in the statistical sense. As a result the binary interaction 
of the phase oscillators looks in the usual way and relates hierarchically to reduced probabilities 
with a larger number n > 2. However, the interaction appears as averaged over the frequency-degree 
distribution (cf. Eq. ([^ ). 
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In the thermodynamic limit N ^ oo the ratio ^iki/{N — 1) becomes the average degree. 
Moreover, the correlation of phases between any two oscillators can be discarded [sj within the 
scope of the mean- field approximation, Eq. @: 

P2 {(t>i,(f>2,t;uJi,uj2,ki,k2) = Pi (0i,t;wi,fci)pi {(f>2,t;uj2,k2) ■ (16) 

Such a decoupling has also been used in [S^, [l^, |3^ [s^ and it is an omnipresent approximation in 
the theory of phase transitions [s^ . A completely rigorous discussion about the general conditions 
has been proposed in ^3§\ by the help of path integral methods. We point out that this effective 
decoupling is fully equivalent to the former description by the help of the stochastic differential 
equation for the single phase oscillator imbedded in the mean field (cf. Eq. 
Let us return to the conditional probability density 

/ , ,, , X Pi ((/>i,t;t^i,fci) , . 

pi{(f>i,t\ui,ki) ^ — ^^^^ — , (17) 

by which we get eventually a nonlinear Fokker- Planck equation for the one-oscillator probability 
density pi: 

dt " ^ b('/'i' Pi {4>iA'^iM)\ + O . (18) 

Therein the mean increment of the phase per unit time reads 

v{(i>\,t) = + rx/ci sin(9 — (/)i) . (19) 
It depends on the density pi(0i, /ci) via the order parameter 

re'"" ^JT7 dcf>2 diU2 / dfcs e^*^ p, {(p2:t\uJ2M) ^2 P(w2,fe) . (20) 

Here >c :— k/N and m > 1 is the lowest possible degree in our network. Eqs. (|18 p - (|^D|) constitute 
again an effective one-oscillator description as Eqs. (I7|), ([5]). Both descriptions are equivalent for 
all possible pairs uoi and fc^, the difference is that in this section the equations are derived for the 
thermodynamic limit. We remark that our expressions are reminiscent of the equations in the work 
of Ichinomiya j22| . who considers uncorrelated random networks without noise. We emphasize that 
the above averaging requires a coarse grained statistical approach to all the nodes. Since each 
of the A'' oscillators is coupled to the mean field with a characteristic strength, we refer to the 
"weighted mean field", which makes sense if the initial conditions of all oscillators are identically 
distributed. We also note that the conditional probability density pi (0i, i|a;i, /ci) in Eq. ([T7]) is 
averaged over the initial conditions. If one would formulate the theory for the transition probability 
pi (01, iji/)'', wi, fci) the nonlinear FPE (|18p would be the same. A difference would occur in the 
definition of the order parameter, where one has to remember that the definition (|20p includes the 
average over the initial conditions as it was performed in Eq. ([l3 



V. THE CRITICAL COUPLING STRENGTH 

To underline the statistical identity of the oscillators in the one-oscillator description, we omit 
the indices and proceed with variables (j),uj,k. Taking into account the normalization condition 

27r 

/3(0,i|w,fc)d(/) = 1 , y oj,k,t , (21) 
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the completely asynchronous stationary state is given by 

po{(f>\uj,k) = —— = const. , W u},k,t (22) 
zn 

with order parameter r — rQ — and undefined phase. It is easy to see that pa is at least one of 
the possible asymptotic solutions of the nonlinear FPE (|20p . 

In the following, we will perform a linear stability analysis of this incoherent solution. In doing 
so, we refer to the work of Strogatz and MiroUo [2^, who also study a nonlinear FPE in order to 
investigate the linear stability of the incoherent state for an unweighted all-to-all connectivity. The 
critical condition, where the incoherent solution loses its stability, is equivalent to the onset of syn- 
chronization we are searching for. To start with, we consider the evolution of a small perturbation 
of the incoherent state: 

p(0, t\oj, k) = ^+ e5p{(j>, t\w, fc), e < 1. (23) 
The normalization condition (|2ip implies 



27r 



(5p(0,t|w,fc)d(/) = V , (24) 



so Sp{(j>,t\u}, k) is 27r-periodic in (p. Now we substitute Eq. ([^5]) into Eq. 

dSp d 



^' dt dcj) 



^+eSp]v 



The amplitude becomes r{t) — eSr{t), where Sr{t) is given by substituting p by 5p in Eq. (pOll . The 
FPE (|18p linearized in the lowest order in e then reads 



dSp 35 p xk ^ ,^ ^9^<5p 

= — tj 1 " — ' 

dt 90 27r 



+ — cos(e -</>)+ D-^ . (26) 



Since 8p is a real number and 27r-periodic in (cf. Eq. (|24p ). it may be expanded as 

oo 

1=1 

Due to Eq. ([M]) . co(i|a;,fc) vanishes. The first nonvanishing coefficients, which constitute the 
so-called fundamental mode, determine the deviations of the mean field in first order of small e: 

Sre'^ = — dio' dk'cl{t\uj',k')k'P{uj\k') . (28) 

\"'/ J — oo J rn 

In order that the incoherent state becomes unstable and a synchronization process starts, ci{t\uj, k) 
has to grow, so that as a consequence the order parameter r{t) grows as well. 

Multiplication of the last equation by e"**^ and considering only the real parts, yields 

6r cos{e - (f>)^^( f du' [ dk'ci{t\uj',k')k'P {uj',k')\ e"^ + C.C. . (29) 

2(k) \J-oo Jm ) 
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Now we insert Eqs. (P7| - (P^ into Eq. and consider only the coefficients witfi e"^. Afterwards, 
one finds tlie evolution equation for the fundamental mode ci{t\uj, k): 

I2l = _(z^ + ,^)ci+ / dw' / Ak'cxW X)^'P{J X) ■ (30) 



dt ' ' 2{k) 

Since this is a partial differential equation without mixed derivatives, we make a separation ansatz: 

ci(t|w, k) = b{uj, k) e^*, A e C. (31) 

So we obtain 

Xb^ -{D + iuj)b+ — duo' I dk'b{uj',k')k'P{uj',k') . (32) 

2(k) J_oo Jm 

The right-hand side is a linear transformation and it can be shown that we only have to calculate 
its point spectrum in order to obtain the critical coupling strength [2^. 

The second term on the right-hand side of Eq. ([5^ equals the degree k times a constant which 
we call B in the following, so that b{uj, k) is given by 

b{u^,k)^ x/n!- ■ (33) 
X + D + 

Hence, we can set up the following self-consistent equation: 

2(k) J -co Jm X + D + iu' 

This relation holds, if all the oscillators are identical in the statistical sense, see section IIVI We 
remark that Restrepo, Ott and Hunt (2005) derived a similar formula where fluctuations emerge 
due to finite-size effects [IJ. We have defined in Eq. ([3]) the value D. As outlined, it does not 
vanish in the thermodynamic limit. In Eq. (|34p the first integral exists for all Re(A) > —D [3^ . 
The solution _B = is not allowed, because otherwise b{uj, k) and as a consequence ci(t|a;, k) would 
vanish. ci{t\uj,k) = however is not considered as an eigenfunction (it is a trivial solution). So B 
is canceled out and we get 



+ 00 



k 



, dw' / dfc'- :P{i^\k') . (35) 

On the analogy of the proof in [4(3|, one can show that there exists only one solution for A and 
that it has to be a real number. The proof makes use of the assumption that, with regard to the 
w-dependency, P{uj, k) has a single maximum at frequency uj = (due to the rotational symmetry 
in the model the maximum can be located like this) and is symmetric. 
In summary the eigenvalue A is given by 

Note that this equation can only be fulfilled if A > —D, because otherwise the right-hand side would 
be nonpositive. Thus, the eigenvalue A is strictly positive and the incoherent solution cannot be 
linearly stable, if the noise intensity vanishes. 
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At the critical condition A = Ac = the incoherent solution loses its stability: if A > the 
fundamental mode ci of the perturbation 5p of the incoherent solution is linearly unstable. It 
grows exponentially with time ^ e"^* and so does the order parameter r{t) (cf. Eq. (l28p V Hence, 
the critical coupling strength >Cc for the onset of synchronization reads 



dw' 



dfc' i: ., p{io\k') 



D2 



(37) 



We emphasize that this equation is not valid in the noise-free case, where one has to take the limit 
A — > 0+ in Eq. with = 0, by which previous results can be reproduced [l^, [2^ (compare 
Tab. IB, 

In [41] we investigate the effects of correlations between the degrees k and the frequencies u}. 
Here we continue with two separated distributions g{uj) and P{k). It is interesting to see that the 
critical coupling strength can then be written as a product of two functionals. The first one 
maps the degree distribution P{k) to a real number via the first two moments, 



/top [P] 



N 



(38) 



We call this one the "topology functional". The second one maps the frequency distribution g{uj) 
to a real number via an integral that depends on the noise intensity D, 



4-00 



D 



We call this one the "diversity functional". In short: 



(39) 



(40) 



VI. APPLICATION TO A DENSE SMALL- WORLD NETWORK MODEL 



In this section we test the analytical expression for the critical coupling strength as a function 
of diversity, noise and network topology for local ring-like networks, where random shortcuts to 
other nodes will be established. For this purpose, we compare the analytical result with numerical 
simulations of networks with different size N. 



A. Derivation of the Topology Function 

First, we specify our dense small- world network model. According to the definition in (4^ . in 
dense networks of TV nodes the total number of edges scales with iV^. Our dense small- world 
networks are constructed as follows. Each oscillator is coupled to its K next neighbors in both 
directions of the ring network and K is given by 



K 



|(iV-l)J , 0<a< 1. (41) 



Here the intensive variable a gives the fraction of all nodes that are coupled through next neighbor 
connections, a = stands for uncoupled nodes, whereas in case of a = 1 the network is fully 
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FIG. 2: (Color online) Visualization of a dense small-world network model. If the system size A'^ is dupli- 
cated, the average number of connections is duplicated as well. The figure shows one marked node in the 
network. 

connected. We underline that the variance of a vanishes. The floor function [.J gives the greatest 
integer less than or equal to its argument, by which the continuous variable a is mapped to the 
discrete variable K. So far we obtain a regular network and due to self-coupling, the fixed degree in 
such regular networks is 2K + 1 . The corresponding degree distribution is equal to the Kronecker 
symbol 

fiocal(fc) = Sk,2K+l ■ (42) 

In our final simulations we choose a = 0.05. Hence, 2.5 percent of possible edges are local regular 
connections. 

Besides these {2K + 1) ■ N next neighbor couplings we introduce shortcuts: each oscillator is 
coupled to the remaining other A'^ — 2K — 1 ones with a certain probability p (see Fig. [5] for 
illustration). If K equals zero the model generates Erdos-Renyi random networks. The case p = 1 
leads to all-to-all connectivity other than in popular sparse small- world network models [ij, . 
The procedure of adding shortcuts can be seen as a Bernoulli experiment with probability p of 
success [27l - l29| . so that the degree distribution P{k) is given by the following binomial distribution: 

p{k) = (I:2^-i)^'"''"[^-^]''"'' 

Note that the minimum degree k ~ 2K + 1 corresponds to Eq. (|42|) . The first moment equals 

{k) = 2K +1 + {N -2K - l)p w (TV - l)(a +p-ap) + l . (44) 

In the last step we have used Eq. (I4ip by neglecting the floor function as an approximation. For 
large N the latter expression approaches {k) k. N{a + p— ap). 
By simple substitutions one finds the second moment 



{k^) = {kf + m-2K-i){i-p) , 



(45) 
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FIG. 3: (Color online) Degree distribution P{k) of dense small-world networks consisting of 400 nodes 
depicted for different shortcut probabilities p and fractions of regular connections a (inset). For the sake 
of clarity solid lines are chosen. 



where the second item equals the variance of k in the binomial distribution due to the randomness 
of the number of shortcuts. It scales only linearly with the system size, whereas the first item 
grows with N"^. Thus, for large systems the difference between the nonrandom number of local 
connections and the number of shortcuts disappears, because the variability of the shortcuts does 
not count for large N. The second moment approaches (/c)^ and it becomes symmetric in a and p. 
Hence, the topology function (cf. Eq. (|38p l for large N results in 



/top(P, a)lAr»i 

We obtain as limiting cases (compare Fig. HJ 



N 



1 



a + p — ap 



/top(p, a) I 



^1 i + ^P + 0{p^), p«l, 
^ l2-a+(a-l)p + 0((l-p)2), p<l. 



(46) 



(47) 



Due to the symmetry, we obtain qualitatively the same for a ^ 1 and a < 1, but with p and a 
being interchanged. As expected, the topology function tends to unity for p — > 1 or a — >■ 1, because 
in both cases the network becomes fully connected. 

In Fig. [3]the dependence of the topology function on a and p is depicted. Discrepancies between 
numerical calculations and theory in the right panel are due to the fact that the number of coupled 
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FIG. 4: (Color online) Dependency of the topology function on the system size A'': Solid lines are given by 
the theory. The thick line corresponds to the thermodynamic limit and slight lines are calculated according 
to Eqs. (|44p and (|45[) for smaller systems. Markers show results of numerical calculations for indicated 
system sizes. 



next neighbors is, of course, an integer. Rounding off in Eq. (|^T|) leads to noticeable steps in 
/top(p, a) for smaller N. 

We emphasize that by the help of Eqs. ([33]) and ([IS]) . one can find the exact expression for the 
topology function for arbitrary system sizes. Since it is a rather lengthy expression, we skip it in the 
text but we will use it in Fig. HI Moreover, as can be seen in Fig. IH a system size N of the order 
O(IOO) is sufRcient to obtain a dynamical behavior that is comparable with the thermodynamic 
limit. Further, the critical coupling strength >Cc has been derived in the thermodynamic limit (see 
Eq. (|4Up ). For this reason it is only consistent to calculate the topology function in the limit 
CX3 as well. 



B. Comparison of the Critical Coupling Strength; Simulations vs. Theory 



In our simulations, the stochastic differential equations ITJ are integrated up to t = 600 with time 
step h = 0.05 by using the Heun scheme [i^. The periods of the oscillators are T ~ 0(10), so our 
integrations cover 0(10) periods. Moreover, in order to calculate statistical equilibria, we discard 
the data up to i = 200, by which transient effects are safely avoided. The statistical equilibria are 
further calculated as averages over at least 100 different network realizations. We emphasize that 
all the different network configurations do not differ only in the configuration of the connections, 
but the oscillators on the network differ as well: all the natural frequencies and the initial values 
of the phases change from one configuration to another one. 

Hong, Choi and Kim (2002) consider the Kuramoto model without noise on a small- world network, 
which is constructed according to the Watts-Strogatz model [i^, where shortcuts are the result of 
rewiring the edges of the initial regular network with a certain probability p [2l[ . They calculate 
the order parameter averaged over time (• • • ) and network realizations [• • • ] and use the following 
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FIG. 5: (Color online) Numerical determination of critical coupling strength by plotting the order parameter 
s times N'^'^^ for different system sizes in dependence on the coupling strength. According to Eq. H48p the 
intersection of the curves indicates the critical coupling strength. Parameter values here: a — 0.05, p — 0.1, 
D = 0.05, a = 0.02. 



scaling form: 



1 ^ 
N ^ 



= N--F 



(48) 



where F[.] is some scaling function. In particular it is found that j3 and v have the same values as 
in the globally connected network, namely P ~ \ and v k, 2. 

The findings propose a finite-size scaling analysis to calculate the critical coupling strength, 
because at >c ~ Xc the function F[.] is independent of the system size N. By plotting s ■ iV°-^^ as 
a function of x for various network sizes N ^ we can measure the critical coupling strength Xc as a 
well-defined intersection point, see Fig. [S] 

In order to compare the theory, expressed by the critical coupling strength (Eq. (|40p l. with 
simulations on the proposed dense small- world networks, we use as topology functional the approx- 
imation for large systems (cf. Eq. (|4t)P ). 

In the following we further consider a Gaussian frequency distribution (7gauss('^) with vanishing 
mean and standard deviation <j. Then for the diversity functional (j39p the expression 



/div(i^)[5, 



gauss 



= 2\ -a 



1 - $ 



D 

7^ 



(49) 



follows; $(.) is the error function. The formula is equivalent to the mean field expression derived 
in [25I I by means of an eigenvalue analysis. In appendix [B] we summarize the limiting cases of Eq. 
(j49p and we give a comparison with other frequency distributions. 

As can be seen in Figs. |n]and[7]we obtain a satisfying agreement for the critical coupling strength 
>fc, regardless of whether we consider the dependencies on the topology (Fig. ^ or on the diversity 

(Fig. m- 

Especially for p > 0.1, a > 0.1 or cr < 0.5 and for the dependency on the noise intensity D in 
general, we obtain almost a perfect agreement between theory and simulations. For smaller values 
of p or a there is a small discrepancy, but the shape of the curves can be well reproduced. 
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FIG. 7: (Color online) The dependencies of the critical coupling strength on the noise intensity D and 
on the diversity parameter a are depicted. Remaining parameters: a = 0.05, p = 0.01. 



As expected, both weaker connectivity and larger diversity impede synchronization. Compare 
Eq. ([T7)) and Tab. U for a summary of the Umiting cases. In particular, for p — > or a — > our 
results suggest a saturation at the values >tc — ^'""^^-'^^J or = li^^d£}M.^ respectively. 

Apparently, the weighted mean-field approximation tends to overestimate the critical coupling 
strength, which is counterintuitive, because the approximation corresponds to a weighted fully 
connected network and all-to-all connectivity should reduce the critical coupling strength. To see 
this, compare p = 1 or a = 1 in Fig. [6l because both choices stand for all-to-all connectivity. So 
the coupling weights defined in Eq. ^ are able to mimic the original complexity in an overstated 
manner. 

The disagreement for high standard deviations a of the Gaussian frequency distribution seems 
to be analogous to the disagreement in the dependency on the topology, because in our derivation 
of the critical coupling strength, both the frequency distribution g{uj) and the degree distribution 
P{k) are involved in averaging. 
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VII. BEYOND THE DENSE SMALL- WORLD NETWORKS 



Since the mean-field description corresponds to all-to-all connectivity, we expect that the approxi- 
mation loses validity in sparser networks. In order to effectively analyze the limitations, we consider 
separately Erdos-Renyi like random networks and regular networks. The first are constructed by 
assigning an edge probability of 

p,=p-m-\ 0<p,g< 1 (50) 

for any two of the A'" nodes in the network. The only additional requirement is that besides the edge 
probability, each node connects a priori to another randomly chosen one. In this way we guarantee 
that there are no isolated nodes, which are not of interest here, because they are not able to take 
part in the synchronization process and only reduce the effective system size. 

In the regular networks each node is coupled to the K next neighbors in both directions of the 
ring network and K is given by (compare Eq. (HH) 



K 



i[3 + a.(7V-4n 



0<a,q<l. (51) 



As for the random networks, we require that there are no isolated nodes, which is already imple- 
mented in the above choice {K > 1). 

In both cases g is a denseness parameter, for q ~ 1 leads to dense and q = to sparse networks 

0. 

We repeat the former analysis, see section lVl Bl with adjusting the normalization of the coupling, 
i.e., we substitute N — >■ TV in front of the coupling term (cf. Eq. ([1])). Consequently, the topology 
functional becomes 

/top[P] = A^''^. (52) 

In case of the regular networks it appears necessary to increase the system size for sparser networks 
in order to have distinguishable networks (see Eq. (HJ)), e.g. for a = 0.001 and q = 0.75 we 
consider network sizes up to iV = 135000. Again, p = 10"-* and a 10"* are the smallest values 
considered in our numerical simulations. 

In Fig. [5]the results are depicted. Unlike below the lines, we observe a mean-field synchronization 
transition above them. In particular, the markers correspond with the smallest p or a values for 
which we observe a synchronization transition with the critical mean-field exponents P — and 
V = 2 [21]. Those p and a values are denoted by p^f or Qfmf- 

As expected, for sparser networks the mean-field approximation breaks down. Furthermore, it 
can be seen that random shortcuts favor the mean-field behavior, since we find the scaling relations 
Pmf ~ 10-^9 and ami ~ lO"^*. 

We also underline that our approach works, if the first two moments of the degree distribution 
exist. So, generally, scale-free networks that display a power-law decay P{k) ^ k~'^ are excluded 
from our approximation. In [23| it was numerically shown that for 7 < 3 the mean-field approach 
yields significant deviations. It is further analytically derived in [3I] that only for 7 > 5 the order 
parameter fulfills a finite-size scaling with the standard mean-field exponenents, and for 2 < 7 < 3 
the exponents depend on the degree exponent 7. Hence, more homogeneous degree distributions 
that can be characterized by the first two moments, favor the mean-field treatment as it was the 
case in our small- world example, see Fig. |31 
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FIG. 8: (Color online) Smallest p (random network) and a (regular network) values for which a synchro- 
nization transition with the standard mean-field exponents can be observed, are shown as functions of the 
denseness parameter q. We denote these values by pmf and Qmf ■ Dashed lines depict the scaling behaviors. 



VIII. CONCLUSION 



We have investigated noisy Kuramoto oscillators on undirected ring networks, whose complex 
structure is approximated by a weighted fully connected network. The weights are obtained by a 
specific combinatorial consideration of the connectivity, where the original network is considered 
to be random and uncorrelated. We have shown that this procedure leads to a weighted mean-field 
approximation, which enabled us to evaluate analytically the critical coupling strength >tf. that 
marks the onset of synchronization in the network. As a result, we have found that is a product 
of two functionals. The first one is a functional of the degree distribution and therefore depends 
solely on the network topology, while the second one is a functional of the frequency distribution 
and a function of the noise intensity. 

As such, we have provided support for the previous separate consideration of network complexity 
and noise with regard to the Kuramoto model in the literature. 

In order to compare the theory with simulations, we have investigated a dense small- world network 
model. As typical of small- world networks [4^1, we start with a regular network consisting only 
of next neighbor couplings, and then we randomly add shortcuts. In this way one can study the 
impact and the interplay of regular local edges and random shortcuts. We have found that in a 
large small-world network, regular local edges and random shortcuts play almost the same role. 

Our dense small-world network model allows scaling between a locally and a globally connected 
network, which is well-suited for testing the mean-field approximation. In case of all-to-all con- 
nectivity the mean- field description is exact; by decreasing the connectivity we can determine the 
arising discrepancies between the original complexity and our estimate. 

It has turned out that the analytical and the numerical results are consistent to each other in par- 
ticular, finite-size scaling analysis with mean-field exponents always gives a well-defined intersection 
point. Hence, the extended Kuramoto model on dense small-world networks shows a mean-field 
synchronization transition. 

We have completed our work by addressing further network models. There we have shown that 
random connections favor a mean-field approach in the sense that regular networks have to be more 
densely connected in order to yield a mean-field synchronization transition. Highly heterogeneous 
networks in which the mean-field description breaks down, have also been discussed. 
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Future work has to be done with regard to network complexity. For instance whether networks 
of higher dimension show significant quahtative differences. In particular, the problem of time- 
dependent coupling strengths or even time-dependent number of nodes and edges has to be tackled. 
A next step could also be the consideration of time delays in the coupling. Furthermore, network 
costs have to be taken into account, which should allow to evaluate the efficiency of a synchronization 
process. 
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Appendix A: Dense Small- World Networks 

Here we provide numerical evidence for the existence of dense small-world networks, see Fig. |H1 
In order to do this, we have to compare such dense small-world networks with random networks. 
However, the random network used for the comparison must not have any isolated nodes. Otherwise 
path length and clustering coefficient are not (well) defined. We take as an almost random network 
(similar to Watts and Strogatz fill) a network generated with our model, but with a minimum 
amount of coupled next neighbors, so we choose K = 1 or an equivalent a (cf. Eq. (|¥T|) ). The 
comparison makes sense, only if the average degree (fc) of both networks is the same, which requires 
choosing the following shortcut probability p' for the random network: 

, 2(K-1) N-2K-1 , 

P = + P ■ (Al 

^ N-3 N-3 ^ ' 

p is the shortcut probability used for the dense small- world network. 

Having in mind the defining properties of a small-world network: L > -Liandom and C ^ Crandom 
[43j . we observe small-world networks for values of p around 10""^. For smaller values of a, i.e., fewer 
coupled next neighbors, the defining properties are better fulfilled. If a is too large the network 
under consideration is no small-world network anymore. 



Appendix B: Diversity Functional for Different Frequency Distributions 



In what follows, we consider the diversity functional for different frequency distributions: uniform 
distribution .guni('^)j identical oscillators 5idcnt(^) and Lorentzian distribution 5iorcntz(w). We obtain 
the following expressions for the diversity functionals (cf. Eq. ([55)1 ): 



/div(-D)[guni] = 2\/3cr 



arctan 



V D 



/div(i5)[5idc„t] =2Z?, 

/div(-D)[giorcntz] =2(i? + 7). 



(Bl) 



a and 7 are the standard deviation and the scale parameter, respectively, while D gives the 
noise intensity. In Fig. [TO] the dependencies on D, cr or 7 are depicted. We observe a greater 
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FIG. 9: (Color online) Numerically calculated average path lengths L and clustering coefficients C as a 
function of shortcut probability p for a network of 501 nodes. The radius of coupled next neighbors (compare 
Eq. (|4ip ') equals K — 2 (Fig. 9(a) I or A" = 50 (Fig. 9(b) I, respectively. Dashed lines indicate the area, in 
which our model generates a small-world network. 



difference in the dependency on diversity than on noise. In particular, for large values of D, all 



the different diversity functional show the same linear dependency on D (cf. Fig. 10(a)). In case 
of Gaussian, uniformly and identically distributed natural frequencies w, it can be shown that the 
diversity functional even approaches the same line for Z? — > cxd. 

However, the above observation is not surprising, because the noise acts on the natural fre- 
quencies. So if the noise intensity D is very high, it makes almost no difference how the natural 
frequencies are distributed, because the diversity of the oscillators mainly comes from the random 
fluctuations induced by the noise. 

Instead, for D ~ const, and increasing diversity parameter a or 7, the different nature of the 
various frequency distributions manifests more and more in a different synchronizability (cf. Fig. 
|10(b)| . Especially the comparison between the uniform and the Lorentzian distribution for the 
same a and 7 values is interesting. In terms of synchronizability we observe that for lower diversity 
the uniform frequency distribution is more favorable than the Lorentzian distribution, while for 
higher diversity the Lorentzian distribution is more favorable. 

Of course, vanishing diversity, i.e., cr — > or 7 — ?> 0, results in identical oscillators, so that the 
diversity functional approaches the same value for every frequency distribution. In Tab. U we 
summarize these results. 



LIMITING CASE 


D < 1 


a < 1 


D> 1 


> 1 


/div(-D)[Sgauss] 

,fd:v(i3)[<?u„i] 
/div(-D)[ffidcnt] 
/div(-D)[fflorcntz] 


2^/^a+^D + 0{D'^) 

2D 
2a + 2D 


2D + %a^ + 0{a*) 

2D + %a^ + 0{a*) 
2D 
2D + 2a 


2Z? + 0(i) 

2D + 0{j^) 

2D 
2D + 0(1) 


2y|a + 0(l) 

^a + 0{l) 
2D 
2a + 0{l) 



TABLE I: The limiting cases of the diversity functional for different frequency distributions. Here a = X. 
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FIG. 10: (Color online) The diversity functional for different frequency distributions as a function of noise 
intensity D, standard deviation a or scale parameter 7. 
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